Finite-size scaling of synchronized oscillation on complex networks 
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The onset of synchronization in a system of random frequency oscillators coupled through a 
random network is investigated. Using a mean-field approximation, we characterize sample-to- 
sample fluctuations for networks of finite size, and derive the corresponding scaling properties in 
the critical region. For scale-free networks with the degree distribution P{k) ^ k ' at large k, 
we found that the finite size exponent v takes on the value 5/2 when 7 > 5, the same as in the 
globally coupled Kuramoto model. For highly heterogeneous networks (3 < 7 < 5), £" and the order 
parameter exponent /3 depend on 7. The analytic expressions for these exponents obtained from 
the mean field theory are shown to be in excellent agreement with data from extensive numerical 
simulations. 
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I. INTRODUCTION 

The popularity of complex networks in the descrip- 
tion of interactions among individuals in various biologi- 
cal and social contexts has inspired theoretical studies of 
ordering phenomena on networks in recent years pj] . The 
small- world properties of such networks, as emphasized 
first by Watts and Strogatz [23, suggest that a simple 
mean-field (MF) description of the ordering transition 
is often appropriate [l|, 0] ■ More complicated situations 
may arise as in, e.g., scale-free networks with a low degree 
exponent, where heterogeneity in the network topology 
smears out the transition significantly [1, 0, H, 0, S Hi • In 
general, randomness in network connections, which can 
be considered as a form of quenched disorder, can have a 
profound effect on the ordering process and the ensuing 
scaling behavior. This is a topic in the network research 
which has not been sufficiently explored so far. 

Synchronization of coupled oscillators is a representa- 
tive dynamical problem on complex networks Jj. By 
varying the coupling strength J among the oscillators, 
various dynamical phenomena can be observed, ranging 
from independent oscillators at J = to fully synchro- 
nized state dX J > Jf. The desynchronization threshold 
Jf is an important property in many applications and 
its dependence on the network topology has been inves- 
tigated in great detail [§]. For coupled oscillators with 
a distribution of intrinsic frequencies [l^, [O, [H, [11] , a 
finite fraction of the population can remain entrained 
in frequency even in the presence of a large number of 
run-away oscillators. The entrained cluster of oscillators 
disappears only at a (much) lower coupling strength Jc- 
Near the entrainment threshold Jc, strong fluctuations in 
various static and dynamic properties of the system are 
expected, as in usual critical phenomena. 

The entrainment transition on scale- free networks with 
the degree distribution P{k) ^ k~'^ has been treated ana- 
lytically by several groups [3, [H, . It has been shown 
that, in the infinite size limit, the transition is expected 



at a finite coupling strength for 7 > 3, while the en- 
trained cluster persists at any nonzero coupling strength 
for 7 < 3. The critical exponent (3 describing vanishing 
behavior of the order parameter on the supercritical side 
is shown to be equal to 1/2 for 7 > 5 and 1/(7 — 3) for 
3 < 7 < 5 [H. 

For systems of finite size, the entrainment transition 
becomes blurred and rounded over a range of the cou- 
pling strength. In addition, randomness in the network 
topology, as well as the random choice of oscillator fre- 
quencies, introduces sample-to-sample variations in the 
entrainment threshold. A full description of the finite 
size effects, which requires a detailed characterization of 
dynamic fiuctuations in specific samples, is currently not 
available [13, [III ■ Fortunately, in the case of the globally 
coupled Kuramoto model where the mean-field theory 
works well on the supercritical side, the sample-to-sample 
fluctuations of the order parameter can be characterized 
analytically . Comparison with numerical simulations 
indicates that temporal fluctuations of the order parame- 
ter only play a subdominant role in the broadening of the 
transition region due to finite size [H, [111 . The success of 
this approach suggests a novel procedure to derive finite- 
size scaling (FSS) relations under a MF approximation. 

In the present paper, we extend the above MF treat- 
ment to coupled random frequency oscillators on complex 
networks. Unlike the globally coupled case, the MF equa- 
tions derived here are not expected to be exact due to the 
finite connectivity of individual vertices on the network. 
Nevertheless, as we show below, the FSS exponents ob- 
tained depend only on certain general properties of the 
network. We also present results from extensive simula- 
tions on uncorrelated scale-free networks. The numeri- 
cally determined values of the exponents as a function of 
the degree exponent 7 of the network agree well with the 
analytic predictions. Our study indicates that the FSS at 
the entrainment transition on uncorrelated scale- free net- 
works is also governed by fluctuations in the distribution 
of intrinsic oscillator frequencies, with temporal order pa- 
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rameter fluctuations playing a less dominant role. 

The paper is organized as follows. In Sec. II we intro- 
duce the dynamic model and derive the mean-field equa- 
tions. Section III contains a treatment of the mean- field 
equations in the neighborhood of the entrainment tran- 
sition. A finite-size scaling form for the order parameter 
is obtained. Results from numerical simulations of the 
model are presented in Sec. IV and compared with the 
analytic predictions. Section V contains a brief summary 
of our results. 



II. THE MODEL AND MEAN-FIELD 
EQUATIONS 

A. The model 

An undirected network of TV vertices is defined by the 
adjacency matrix {cty }, where Uij = 1 if two vertices i 
and j are connected and otherwise. The degree of a 
vertex i is the number of vertices connected to i, denoted 
by ki = J2j (^ij ■ To each vertex i we associate an os- 
cillator whose dynamics is described by the equation of 
motion for the phase 0^, 

N 

4>i = uJi - j'^Uij sm{(j)i - (pj), (1) 
i=i 

where uji is the intrinsic frequency of the oscillator. The 
second term on the right-hand side of Eq. ^ denotes 
coupling to neighboring oscillators on the network with 
a positive strength (J > 0). 

In this paper, the oscillator frequencies uji in a given 
sample are drawn independently from a distribution g{uj) 
which is assumed to be a smooth function and symmetric 
about its maximum at w = 0. In addition, we shall limit 
ourselves to random networks with no degree-degree cor- 
relation among neighboring vertices. An algorithm that 
generates such a network with no self linking nor multi- 
ple links between vertices is discussed in Ref. [23|. The 
numerical results presented in Sec. IV are for scale-free 
networks generated using the static model of Ref. ^] . 

B. The mean- field approximation 

At sufficiently strong coupling, the system described 
by Eq. ([1]) exhibits a synchronization phenomenon where 
a finite fraction of oscillators in the system become en- 
trained in frequency. Transition to the random state at 
a critical coupling strength Jc has been considered at the 
mean- field level by several authors [3,[l^,[lB|- To under- 
stand the general idea behind such an approach, let us 
first introduce a set of instantaneous local fields defined 
by 

F,e^^' =^a,,e**^ (2) 
j 



where Hi > denotes the amplitude and 0i the phase, 
respectively. Eq. ([T]) can now be written in a more sug- 
gestive form, 

(pi = LOi - HiJ sm{(l)i - 0i). (3) 

As usual, the mean-field approximation decouples the 
set of N dynamical equations by replacing Hi and 9i with 
suitable time-averaged quantities which are then deter- 
mined in a self-consistent manner. Unlike the globally 
coupled Kuramoto model, such a mean-field treatment is 
not exact due to the finite connectivity of individual oscil- 
lators on the network, so dynamic fluctuations are not av- 
eraged out even for an infinite system. However, provided 
the network is sufficiently well connected (as compared 
to fragmentation into local "communities" ) , and because 
of its small-world properties, once a cluster of entrained 
oscillators is formed, it will affect the dynamics of all os- 
cillators in the system through a global "ordering field" 
H ex-p{i9). The precise effect of the global field on indi- 
vidual oscillators is subject to renormalization by the lo- 
cal community of a given oscillator. For a well-connected 
network, it is reasonable to expect that renormalization 
of oscillator's response function does not change qualita- 
tively results of the mean- field theory where such effects 
are ignored. This is to be confirmed by numerical inves- 
tigations. 

With the above caveat, let us proceed to the derivation 
of the mean-field equations, with particular emphasis on 
the finite size effect. The key approximation we introduce 
is to replace each phase factor e"^J in the sum of Eq. ^ 
by the global ordering field H exp{i9) acting through the 
edge connecting vertices i and j. Equation ^ now be- 
comes 

(j)i^uj^- k^HJsm{(j)i - 9). (4) 

Consider a steady-state situation with a constant H 
and a linearly advancing 9 — fit, where i7 is the phase 
velocity of the entrained cluster. From Eq. (|4]), the oscil- 
lator at vertex i is entrained with 

cj),^9 + sm~^[{uj^^n)/{hHJ)] (5) 

if l^j — f2| < kiHJ, and detrained otherwise. In the latter 
case, the time averaged value of e*'^"^'"^^ is given by 

Here and elsewhere the overline bar denotes time average. 

The self-consistent equations for H and are obtained 
by setting H equal to the average of e'^^'^'~^'> over all edges 
of the network. Since each vertex i contributes ki times 
to the average, we may write 
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Here (k) is the average degree of a vertex in the network. 
Substituting Eqs. © and ^ into Eq. ([7]), and separating 
out real and imaginary parts, we obtain 



N 



(8) 



N 



= Y^^cj, -n) i-^i- fr^em - 1) , (9) 



where fi = {ui — fl)/{kiHJ), and 8(2;) is the Heaviside 
step function which takes the value 1 for x > and 
otherwise. 

Equations ([5]) and © are our mean-field equations for 
a given network of N oscillators with a particular set of 
intrinsic frequencies {coi}. They are invariant under a 
uniform shift of all oscillator frequencies (with a corre- 
sponding change in Q). 

At any value of J, Eqs. ([S]) and © admit a trivial 
solution H = (with arbitrary fi). In the next section, 
we shall analyze nontrivial solutions that appear above 
the entrainment threshold Jc, with particular attention 
on sample-to-sample variations when N is finite. 



III. SOLUTION OF THE MEAN-FIELD 
EQUATIONS 

From Eqs. ([5]) and (O it is obvious that the solution 
for H depends not only on the coupling strength J but 
also on the particular choice of the intrinsic frequencies 
{wi} as well as the particular network topology in a given 
sample. For sufficiently large N, it is possible to give 
a statistical description of the sample-to-sample fluctu- 
ations. As we show below, the analysis also yields the 
finite-size scaling of the critical properties in the neigh- 
borhood of the entrainment transition. By symmetry we 
shall seek a solution at 17 = and ignore weak finite-size 
corrections which do not alter our main conclusions. 



A. Self-averaging in the infinite size limit 

To proceed, let us write Eq. ([H]) in a symbolic form, 

H^^iH), (10) 

where 



(11) 



Terms in the sum can be grouped according to their de- 
gree ki- When the network size N 00, the number of 
terms in each group at a given k grows linearly with N, 
and hence the self-averaging over the distribution g(uj) is 



expected. This consideration leads to the result 

1 



«'(i?) = lim *(i?) = (*(iJ)) 



,^.Y' P(k)ku(kHJ). 

(12) 

Here P{k) is the degree distribution of vertices on the 
network, and 



u{x) 



dujg{uj) \/l — up- jx^ 



(13) 



is a monotonically increasing function of x which ap- 
proaches 1 as a; ^ 00. For small x, u(x) ~ ^g(0)a; -I- 
§5"(0)2;^ 

We now consider the behavior of ^'(i?) at small i7, 
assuming (fc^) = J2k^^^(^) finite (i.e., P{k) falls 

off faster than k^^ at large fc). To facilitate the analysis, 
we write 



u{x) = ax — u{x)^ 



(14) 



where a — ^g{0) and u{x) ~ —j^g"{0)x^ for small x. 
Substituting Eq. (HH) into Eq. we obtain 

^{H)=a^JH-i>iH), (15) 



where 



ii{H) ^4xYl P{k)ku{kHJ). (16) 



The functional form of '^{H) at small H depends on the 
tail of the degree distribution P{k). If P{k) falls faster 
than k~^ at large fc, we may use the small x expansion 
of u{x) to obtain 



*(i^):.co^(Ji^)^ 



(17) 



where cq — —{'!i/lQ)g"{Q) is a positive constant. On the 
other hand, if P{k) ~ Ak^"' at large k with 3 < 7 < 5, 
(fc^) diverges and the above expansion becomes invalid. 
Instead, contributions to the sum in Eq. (jl6p come 
mainly from vertices with k ^ (JH)^^. Since the frac- 
tion of vertices in this degree range is proportional to 
k^"'^^, and each contributing an amount k to the sum, 
we estimate ^{H) ~ {JH)'^~^ in this case. More pre- 
cisely, replacing the sum over k by an integral, we obtain 

^(H) ~ — / dkP(k)ku(kHJ) 
(k) Jo 



1 r°° 

— / dkAk-'-'+^mnj) 

k) Jo 



(k) 
cii-JH) 



7-2 



(18) 



Here Ci = A{k) ^ dxx '^^^u{x) is another positive 
constant. 
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B. Sample-to-sample variations in a finite networlc 

To work out the statistics of the sample-to-sample fluc- 
tuation S^{H) = ^{H) - ^'(7?), we note that ^{H) can 
be viewed as the mean value of the random variable, 



over N realizations. From the central limit theorem, we 
expect 6^{H) to satisfy the Gaussian distribution with 
zero mean and a variance given by 

mr) = ^{{v')-{vr)^^- (20) 

As {t]) = "^{H) ^ H sd small H, we only need to focus 
on 



(r;2)=^P(fc) 



fc2 



kHJ 



dLog[ijj) 



kHJ 



" KkHJJ 



(21) 



close to the transition. If P{k) falls faster than fc""* at 
large fc, the above expression can be easily evaluated at 
small 7J to give 



D{H) ^ (ry^) 



2^ 4(fc^) 



g{0)JH. 



(22) 



On the other hand, if P{k) ~ Ak with 3 < 7 < 4, we 
may replace the sum over k by an integral which yields 



D{H) ~ di{JH) 



7-3 



(23) 



where di — A{k)^^ duv?^^ J^^ dujg{uj){l — uj^/u'^) is 
another positive constant. 

Summarizing the above results, in the neighborhood 
of the entrainment transition, we may write the self- 
consistent equation (fTU|) for H in the form, 

H = {J/Jc)H - c{JHf + d{JHY/^N-^''^i, (24) 

where c and d are positive constants, and Jc — 
2{k) / {■Kg{Q){k'^)) (mean-field value) is the critical cou- 
pling strength at the transition. The term ^ is a Gaus- 
sian random variable with zero mean and unit variance 
that represents the combined effect of random frequency 
and network realizations in a given sample. The expo- 
nents p = 3 and q = 1 when the degree distribution P{k) 
decays sufficiently fast at large k. For power-law distri- 
butions P{k) ~ p switches to the value 7 — 2 for 
7 < 5. Similarly, q switches to the value 7 — 3 for 7 < 4. 



C. Finite-size scaling 



The network size N enters Eq. ([24| through the last 
, where ^ varies from sample to sample. At J — Jc, a 
positive ^ (50% of the samples) yields a solution H ^ 



7V-i/(2p-<?)^ quite a bit larger than the value 7V~^/^ due 
to dynamic fluctuations on the detrained side [1^. In 
addition, one can easily verify that the term ^ affects 
the solution significantly when J is within a distance of 
order 7V~'^p~^^/*^^p^'^) from Jc- These properties are sum- 
marized in the scaling solution to Eq. 

(25) 



where the exponents (3 ~l/{p—\) and D ~ {2p — q)/{p — 
1). Unlike the usual FSS expression for pure systems, the 
scaling function f{x) varies from sample to sample. 

Using the above values for p and q, we obtain the fol- 
lowing results for the critical exponents in case of uncor- 
related scale-free networks. 



7—3 ' 7—3 
1 7-1 

7-3' 7-3' 



7 > 5; 
4 < 7 < 5; (26) 
3 < 7 < 4. 



In simulation studies, the order parameter 



A 



N ^ 



-B) 



(27) 



is usually measured. From the solution to the mean-field 
equation ([?]), we obtain 



1 ^ 



k,HJ 



6 1 



k,HJ 



(28) 



Following the same procedure as above in the calculation 
of ^'(i?), we obtain at small iJ, 



A ~ a{k)JH + d2{JHfl'^N-^''^l, 



(29) 



where ^2 = [(4/3)5(0) (fc)]^/^ and ^ is another Gaussian 
random variable. The finite-size term here is negligible 
even at the transition. Hence the scaling behavior of A 
is the same as that of H . 

We note that the expressions for the FSS exponent v 
in Eq. (pS)) differ from those obtained by Lee [l^ based 
on a cluster analysis but without considering sample-to- 
sample fluctuations. In fact, v obtained here is always 
larger than that in (Tsj for any 7. This observation sug- 
gests that the broadening of the transition region by the 
random nature of oscillator frequencies and network con- 
nectivity is more significant than other effects such as 
those considered in Ref. ^15i] . 



IV. NUMERICAL RESULTS 

To test the validity of our MF analysis, we have per- 
formed extensive numerical simulations on the system 
governed by Eq. Scale- free networks used in the 



simulation are generated following Ref. [24I , up to a sys- 
tem size A'' = 12800. The intrinsic frequencies of oscilla- 
tors are drawn independently from the Gaussian distribu- 
tion g{u}) = (27r(T)~"'^/^ exp(— ti;^/2cr^) with unit variance 
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FIG. 1: (Color online) Critical decay of A at the critical point 
Jc as a function of A'^ for 7 — 7.07, 4.48, 3.75, where Jc = 
1.86(2), 1.57(4), 1.17(3), respectively. Lines are drawn with 
slopes of ~l3/h' given by the theoretical prediction, Eq. (|26|) . 



(cr^ = 1). The Heun's method [131 j with a discrete time 
step 6t — 0.01, is used to integrate numerically Eq. ([1]). 
Typically, the motion is followed for A^t = 4 x 10'' time 
steps, with the initial condition 0^ = for all i. Data 
from the first Nt/2 steps are discarded in measuring the 
time average of various quantities of interest. For each 
network size, approximately 10'* independent runs are 
performed, with different realizations of the intrinsic fre- 
quencies as well as network connectivity to obtain sample 
averages. 

To characterize the entrainment transition numeri- 
cally, we have focused on the order parameter defined 
by Eq. (P7)) . or more precisely, A = (At), where A; — 
jj X^jLi exp(z(^j ) I at time t and (■ • ■ ) denotes the sample 
average for a given size. The critical coupling strength 
Jc is estimated from the crossing point in the plot AiV*^ 
versus J, varying the exponent e. This value for Jc is 
checked against crossing of the Binder cumulants at dif- 
ferent size defined by = 1 - i(Af/Af^) [H, and 
compared with the behavior of the dynamic susceptibil- 
ity X — N{A^ — At ). At the critical point J = Jc, the 
order parameter A exhibits a power law: A ~ N~^^'^, 
which provides an alternative way for checking the criti- 
cal value as well as determining the exponent ratio 
With the values of P/D so obtained, we estimate the value 
of the ESS exponent D from the scaling plots of AN^^'^ 
against (J — Jc)N^/'^ for a broad range of systems sizes 
N, adjusting D to achieve the best data collapse in the 
critical region. 

Eigure [1] shows the critical decay of the order param- 
eter A at Jc as a function of N on log-log scale. The 
critical values of J are given by Jc = 1.86(2), 1.57(4), 
and 1.17(3) for 7 = 7.07,4.48, and 3.75, respectively. 
Each set of data points at a given 7 fall nicely on a 
straight line corresponding to the anticipated power law 
behavior A(Jc) ^ N~^^'^, and the slope shows very good 
agreement with the predicted value for (3/9 according to 
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FIG. 2: (Color online) Scaling plot of A for iV = 
800, 1600, ■ ■ ■ , 12800 with (a) 7 = 7.07, using 13/u = 1/5 and 
u = 5/2, (b) 7 = 4.48, using f3/u = 0.25 and v = 2.68, and 
(c) 7 = 3.75, using P/v — 0.36 and v — 3.67. 



Eq. 1^. 

Eigure [2] shows the scaling plots of the order parameter 
A for various system sizes at three different values of 7. 
In each case, the theoretically predicted values given by 
Eq. ([26]) are used to scale the horizontal and vertical axes. 
The data collapse is nearly perfect for 7 — 7.07 and 4.48, 
and satisfactory for 7 = 3.75. 

Einally, we present in Eig.[3]results for /3/P and P deter- 
mined following the procedure described above at various 
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FIG. 3; (Color online) P/u (box) and l/u (triangle) are plot- 
ted as a function of the degree exponent 7, showing a good 
consistency with the theoretical prediction represented by the 
full line {P/v) and the dotted one (1/j^), respectively. 

values of 7. Error bars are obtained from uncertainties 
in the procedure. The two lines correspond to our MF 
predictions. The agreement is generally good. 

V. CONCLUSIONS AND DISCUSSION 

In this paper, we investigated the effect of fluctuations 
in the frequency distribution on the entrainment transi- 
tion of a system of coupled random frequency oscilla- 
tors. Self-consistent equations for the global ordering 
field are derived for any given oscillator population on 



a complex network. Statistical properties of the solution 
to these equations in an ensemble of such networks of 
oscillators are determined. The analysis enables us to 
derive a finite-size scaling expression for the entrainment 
order parameter. Comparison with numerical integration 
of dynamical equations on scale-free networks shows that 
the mean-field description correctly captures the finite- 
size scaling behavior exhibited by the sample-averaged 
order parameter A. The values of exponents (3 and v 
which best describe the numerical data show excellent 
agreement with the predicted ones from the mean-field 
theory. 

For scale- free networks with a degree exponent 7 > 5, 
our result for the exponents /3 = 1/2 and v = 5/2 is 
the same as the globally coupled Kuramoto model. Ran- 
domness in the network connection does not appear to 
affect these values. On this ground we expect the re- 
sult to apply to randomly connected networks with a 
bounded degree, to the Erdos and Renyi network [2^, as 
well as to the small- world networks of Watts and Stro- 
gatz [3]. However, the finite connectivity of oscillators 
on the network implies that the dynamic fluctuations are 
not averaged out and can at least renormalize the en- 
trainment threshold Jc given by the mean-field theory. In 
this regard, analytic derivation of the sample-dependent, 
self-consistent equation ([24]) beyond the mean-field ap- 
proximation would be desirable. 
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